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Clusters generated by the product-rule growth model of Achlioptas, D'Souza, and Spencer on 
a two-dimensional square lattice are shown to obey qualitatively different scaling behavior than 
standard (random growth) percolation. The threshold with unrestricted bond placement (allowing 
loops) is found precisely using several different criteria based upon both moments and wrapping 
probabilities, yielding p c — 0.526565 ± 0.000005, consistent with the recent result of Radicchi and 
Fortunato. The correlation-length exponent v is found to be close to 1. The qualitative difference 
from regular percolation is shown dramatically in the behavior of the percolation probability P^ 
(size of largest cluster), the susceptibility, and of the second moment of finite clusters, where discon- 
tinuities appears at the threshold. The critical cluster-size distribution does not follow a consistent 
power-law for the range of system sizes we study (L < 8192) but may approach a power-law with 
r > 2 for larger L. 

PACS numbers: 64.60.ah, 64.60.De, 05.50.+q 



I. INTRODUCTION 

Recently, there has been a great deal of interest in a 
model of "explosive growth" of percolation clusters by the 
so-called Achlioptas process pQ, in which two randomly 
chosen unoccupied bonds in a system are examined, and 
the bond that minimizes the product of the size of the two 
clusters to which it is attached becomes the next one to 
occupied. This procedure, called the product rule (PR) 
PQ , was originally studied on Erdos-Renyi random graphs 
PQ , then on two-dimensional square lattices [2] and scale- 
free networks [2111]. Other recent papers on explosive and 
biased percolation include [5Ul6| . Interest in this process 
derives from its unusual explosive behavior, suggesting 
a first-order transition, with apparent discontinuities in 
several quantities. Many of its properties have yet to be 
discovered. 

In this paper we examine in more detail the PR model 
on the regular square lattice, especially in regards to how 
it differs from random growth (RG), in which bonds are 
added one at a time, and which corresponds to standard 
percolation. 

Some preliminary results were given in [2 , where the 
width of the distribution A/N was investigated. As in pQ, 
A/N was defined as the difference in times in which the 
maximum cluster size s max goes from to 0.5N, where 
N is the number of sites. (Here time is identical to the 
number of bonds added.) For the Erdos-Renyi graph, the 
authors of pQ found that A / N -» as N -> oo for the PR 
model, while A/N — > const, for the RG model, showing 
that the two transitions are qualitatively different. In 
[2] it was found that for the square lattice, on the other 
hand, A/N -> for both the PR and RG models, but 
with different powers in A^. It however turns out that if 
a larger (and in the case of the square boundary, more 
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appropriate) criterion for the upper end of the gap A 
were used, say s max = Q.7N, then indeed one would find 
that A/N — > const, as iV -> oo for the RG model and 
still to zero for the PR model. So with this criterion, 
the two models are qualitatively different on the square 
lattice just as for the Erdos-Renyi graph. (Even with a 
criterion of s max = 0.5N, one should have A/N — > const, 
for the RG model on the square lattice, but one would 
have to go to a very large system to see it.) 

Recently, Radicchi and Fortunato have also studied 
both the PR and RG models on the square lattice, and 
analyzed their behavior in the context of standard two- 
parameter scaling |14j . However, as they mention, it is 
unclear in what sense this scaling can be applied to the 
PR problem, considering that several of the quantities 
show discontinuities. In this paper, we consider the be- 
havior of wrapping probabilities as well as quantities re- 
lated to the size distribution such as moments. The for- 
mer refers to having a cluster that connects around the 
toroidal boundaries of the periodic system (the torus), 
and is the analog of crossing probabilities for open sys- 
tems. In order to study scaling behavior precisely it is 
also necessary to know the transition point precisely, and 
we determine it here using a variety of methods. While 
the convergence of various estimates in ordinary perco- 
lation is well-known p~7| [18] , that is not the case for the 
PR model, so the convergence behavior is also studied. 



II. PROCEDURE 

Actually, there is a subtle but significant difference in 
the treatment of the PR process considered previously 
by the present author [2] and that by Radicchi and For- 
tunato 14J. In [2], it was assumed that bonds could only 
be added between different clusters. This assumption by- 
passed the question of how to assign weights when a bond 
connects sites that are part of the same cluster, and for 
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FIG. 1. (Color online) (s max )/iV (or Poo) vs. p as a function of system size. Here, as in all of the figures, we show curves for 
L — 128 (red), 256 (orange), 512 (green), 1024 (blue), and 2048 (violet) — in general, from more gradual [L = 128) to sharper 
(L = 2048). Vertical dashed lines show the transition points, p c — 0.5 (RG) and p c = 0.526565 (PR). Plots for RG are always 
shown on the left, and those for PR are shown on the right. The scaling behavior of (s max )/JV at p c is shown in Fig. [2] 
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FIG. 2. (Color online) ln((s max (p c ))/Af) (upper points) and ln(M2(p c ) / N) (lower points) as a function of InL, where p c = 0.5 
(RG) and p c — 0.526565 (PR). Linear fits to the points are shown on the plots, where y represents the the abscissa value and 
x represents InL. In these plots, we have also included data for runs at L = 64. 



the RG model corresponds to "loop-less" percolation pre- 
viously considered by Manna and Subramanian |19j . On 
the other hand, in [U the authors considered that bonds 
can be placed anywhere, including within the same clus- 
ter. Thus, there is a difference in the scaling of the time, 
but also a difference in the weights with which new bonds 
are added, so these two processes are not equivalent. 

In this paper, we follow the unrestricted bond place- 
ment convention used by Radicchi and Fortunato in [14] . 
When an internal bond is selected, we use for its weight 
the square of the size of the cluster it is part of. We 
characterize the size of a cluster (or component) by the 
number of sites it contains. 

To carry out these simulations, we used the algorithm 
of Newman and Ziff [2U1 [H] in which clusters are repre- 
sented as a tree and a union-find algorithm (modified for 
the PR) is used to join clusters together. A randomly 
ordered list of all bonds is made initially, and bonds are 
taken off that list in pairs. The bond that is not selected 
according to the PR is put back on the list randomly by 
switching with a randomly chosen member remaining on 
the list. We also considered the less efficient procedure of 



not using a bond list but just randomly selecting bonds 
on the lattice, skipping over those that were already cho- 
sen until two free ones were found. Both methods led to 
identical results. 

To determine cluster wrapping, we assigned extra vari- 
ables "xcoor" and "ycoor" to each lattice site. These 
quantities are the x— and y— coordinates of that site 
with respect to the first site of the cluster, without ad- 
justing for the periodic boundary conditions. Wrapping 
is then indicated when an intra-cluster bond leads to a 
difference in a coordinate by a multiple of the lattice 
width or height [22] . 

The algorithm of [20] allows one to find the various 
quantities for all values of p in one simulation. We did 
not carry out the convolution step with a binomial dis- 
tribution to get the grand canonical (fixed probability) 
rather than canonical (fixed number) results, as the dif- 
ferences between the two ensembles for the systems we 
studied are small. We everywhere consider square lat- 
tices and square boundaries, with N = L x L sites and 
periodic b.c. Many runs were made to get good statistics, 
ranging from 1000 000 runs for L = 128 to 150 000 runs 
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for L — 2048. In general, the number of runs was suffi- 
cient so that the errors are smaller than the symbols or 
width of the lines or symbols we used to plot the results. 
We also considered runs for L = 8192 for measuring the 
size distribution. 



III. RESULTS 

The results of this work are shown in a series of pairs 
of figures, with results for the RG model placed on the 
left-hand side, and those for the PR model placed on the 
right-hand side. 



A. The maximum cluster size 

In Fig. [I] we show the average of the maximum cluster 
size scaled by the number of sites, (s max )/N, as a func- 
tion of p, for the different system sizes. This quantity 
can also be identified with the usual order parameter, 
the percolation probability , if one considers that the 
largest cluster is effectively the "infinite" one. The PR 
model (right panel) clearly shows qualitatively different 
behavior than the RG model, with crossing curves in the 
PR case. 

The behavior of (s max ) /N at p c is shown in Fig. [5] 
using p c = 1/2 for the RG case and p c = 0.526565 (de- 
termined below) for the PR case. For the RG case, the 
slope (—0.1062) agrees within errors with the scaling pre- 
dictions of -P/v = -5/48 ps -0.104167. The points for 
the PR model are also fit well by a straight line on the 
ln-ln plot, suggesting scaling for this quantity, with slope 
—fi/v = —0.0589, which is clearly different from the RG 
model. Based upon the variation with size, we estimate 
the error to this result to be ±0.01. This value of fi/v is 
consistent with the value f3/v = 0.07(3) (within the error 
bars ±0.03) reported in [13] . 



B. Moments and susceptibility 

Fig. [3] shows the behavior of the second moment 
M%{p) — J2 S s2 "s = (1/^0 J2i s h where S{ is the mass of 
the i-th particle, scaled by N. Again, the PR model 
shows curve-crossing with a possible accumulation or 
crossing point. The scaling behavior at p c is shown in 
Fig. [2] The slope for the RG model —0.210 is consistent 
with the prediction 7/1/ - 2 = -5/24 pa -0.208333. The 
PR data also appears to obey power-law behavior, with 
a slope j/v — 2 pa —0.10 implying 7/1/ pa 1.90, with an 
estimated error of 0.01. This is somewhat higher than 
the value j/u = 1.7(1) reported in [T4"] . 

By scaling and hyperscaling in 2d, one would expect 
that the slopes of the two curves in Fig. [2] should differ 
by a factor of two: 7/z/ — 2 = — 2/3/za This is seen to 
hold well for the RG data, but not so well for the PR 
case. Further analysis of the data shows that the value 



j/v — 2 ss —0.10 seems to be independent of L, but —(3/v 
appears to be increasing as L increases, and may possibly 
approach the value —0.05 (as L — > 00) implied by this 
scaling. However, studies on larger systems are needed 
to confirm this. 

In Fig. |4j we show the behavior of the scaled second 
moment minus the largest cluster, that is: 
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According to scaling arguments, this function should go 
through a maximum at a value p — p max where the 
(Pmax ^ PcjL 1 ' 11 is a certain constant, at which point, 
Mz(p max )/ N should scale as L 1 ^^ 2 . We verified that the 
peaks for RG in Fig. [4] obey this behavior with standard 
exponents (not shown). However, for the PR model, the 
curves of M%/N very closely pivot around the crossing 
point at p c pa 0.52654, which is also close to the inflec- 
tion points of the curves. This suggests that as L — > 00, 
M^Pcj/N is non-zero, which would imply that j/v = 2, 
in conflict with what we found (7/1/ ps 1.90) from the 
behavior of Mi (p c ) . This behavior is another indication 
of the unusual nature of the PR transition. 

In Fig. [5] we show the behavior of the susceptibility x, 
defined by 



X — \/( s max) ( s max) 2 , 



(2) 



which characterizes the fluctuations in the size of the 
largest cluster. It can also be found from previous quan- 
tities via x = [M-2 - M2 - (w) 2 /^) 1 ' 2 - Tne P eaks 
of x{Pc)/N in the RG model decay to zero as L~ - 22 , 
consistent with the scaling prediction L 1 ^^ 2 = L~ 5 / 24 . 
However, the peaks in the PR model apparently increase 
to a constant value x/iV pa 0.264, again consistent with 
7/1/ = 2. Also, the locations of p at the peaks for the PR 
model approach p c pa 0.526575 as L^ 1 , again implying 
v = 1. Below we will find that several other quantities 



also satisfy inverse-size scaling (Fig. 10). 

In Fig. § we show a scaling plot of £/N vs. (p — p c )L 
assuming v — 1, and also 7/1/ = 2, and the fit is seen to 
be good. (Taking v — 1/0.96 yields a much poorer fit.) A 
similar plot for the RG model, with standard percolation 
scaling, is shown for comparison. 



C. Wrapping probabilities 

Next we consider wrapping probabilities. For stan- 
dard percolation these were studied theoretically by Pin- 
son [23] and numerically in [HI [24j [25] . This work has 
also been generalized to the Potts model [3S] [27] . Even 
though the percolating critical cluster is a fractal and of 
zero density in the continuum limit, the wrapping prob- 
ability remains finite and has a value that depends only 
upon the aspect ratio of the system and the type of wrap- 
ping homology. 
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FIG. 3. (Color online) Scaled second moment M2(p)/N as a function of p. 
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FIG. 4. (Color online) Scaled second moment minus largest cluster M2(p)/N, showing a distinct qualitative difference between 
the two models: the curves in the PR model cross at one point (presumably p c ), while those of the RG model do not. Lower 
plots show close-ups around p c . 



In Fig. u\ we show the (only) one-way wrapping prob- 
ability II", defined as the probability at least one clus- 
ter wraps horizontally but not vertically, or wraps ver- 
tically but not horizontally. For the RG model, the 
value of II' 1 ) at p c = 1/2 approaches the predicted value 
0.351642855 . . . [2TJ H3j very rapidly. The curves are ex- 
actly symmetric, because for one-way wrapping there 
must also be a one-way wrapping on the dual lattice, 
which in the square system is identical to the original 
lattice but with bonds occupied with probability 1 — p. 
For the PR model, the curves are not quite symmetric, 
and the value of p at the peaks approaches 0.52658 appar- 



ently as L^ 1 (not shown). The value of seems to be 
dropping to a constant value of about 0.18 as rj L~ 5 , 
although the range of values of L we considered is not 
sufficient to be very certain about this behavior. 

The width (standard deviation) of U^(p), as a func- 
tion of L, is shown in Fig. [8] For RG, the data are consis- 
tent with the theoretical prediction of a straight line with 
slope of — l/v = —0.75. For the PR, the overall slope of 
the points is —0.95 but decreases to —0.96 for large L, 
implying that v w 1/0.96. This is in contrast with the 
value v ~ 1 seen in several other situations. 

In Fig. [9] we show the probability distribution 11^ for 
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FIG. 5. (Color online) Susceptibility \/N as a function of p. Lower plots are close-ups of the behavior near p c . 
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FIG. 6. (Color online) Scaling plots of the susceptibility, assuming v — 4/3, r y/2v = 43/48 (RG), and v = 1, 7/2^ = 1 (PR). 
In both cases, curves L = 128 have the lowest peak, and L = 2048 have the highest peak. 



horizontal wrapping, irrespective of whether wrapping 
occurs in the vertical direction. For both the RG and 
PR models, the curves cross at a single point, within nu- 
merical error. The crossing point of the RG model is 
at p = 0.499995(5), = 0.5210, consistent with Pin- 
son's theoretical result Ii (h) {p c ) = 0.52105829 . . . [2TII215] . 
while that for the PR model is at p = 0.526566(3), 
IlW = 0.5106. Conver gence behavior of the ordinary 
percolation crossing point was studied in [21] . however 
for site percolation and in the grand-canonical (fixed p) 
rather than the canonical (fixed-n) ensemble. We have 
not determined the convergence in this case, but the 
crossing point for the system sizes we considered clearly 



gives a very precise indication of p c . 

The estimates for p c that come from various mea- 
sures are summarized in Fig. [lOj plotted as a function 
of L x l v = L-°- 75 (RG) and L^PR). The upper curves 
show the average of p at which one-way wrapping oc- 
curs — that is, the mean of the distribution shown in 
Fig. [7J ^pIL^^p). The middle curves show the average 
value of p at which horizontal wrapping first occurs; for 
the PR case, this is nearly horizontal, so this quantity is 
very good for estimating p c precisely. The lower curves 
show the average value of p at which either horizontal or 
vertical crossing first occurs. For the RG model, all esti- 
mates extrapolate to a value very close to the expected 
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FIG. 7. (Color online) One-way wrapping probability IT' 1 ' as a function of p. The width of the distribution is plotted in Fig 
[8] and the mean p values are plotted in Fig. |10| 



4-» 



5-5 



rts 



CD " 

c 

o 







RG 








. y = -0.7387x 


- 1.1038 ~ 





6 7 
In L 



S -6 

5 
i 

c 
o 



-7 
-8 
-9 





PR 


<3>^ 












" y = -0.9535X - 1.1178 v 




i i i 





6 In L 7 



8 



FIG. 8. The In of the width of n (1) (p) (shown in Fig. [7} as a function of lnL. The linear fit to the points is shown on the 
plot, where y represents the In of the width and x represents InL. The slope is consistent with 1/v = 3/4 for the RG case, and 
suggests 1/v ~ 0.95 in the PR case. 



value 0.5, and for the PR model the extrapolations are 
consistent with p c — 0.526263(3). 

Note that here, we find a better fit to the data assum- 
ing v = 1 rather than v « 1/0.96 found in the scaling of 
the one-way width, Fig. [Hj However, if we use the latter 
value, we don't find a significant change in the estimated 
value of p c . 



D. Size distribution 

Finally, we consider the behavior of the cluster size 
distribution at criticality. We ran simulations on sys- 
tems of size L = 512, 2048, and 8192 for both the RG 
and PR models, and measured n s , the number of clus- 
ters of size s, at the critical ^points p c . We binned the 
weighted data as P s — Y1J= S s ' n s' f° r s = 1; 2, 4, 8, . . ., 
thus accumulating the number of occupied sites belong- 
ing to clusters in each size range. That is, when growing a 
cluster of size s, we incremented the bin n = (int)(log 2 s) 
by s. (This gives better statistics than the usual method 
of incrementing the bin by 1, which corresponds to just 



counting the number of clusters in each size bin.) For 
a given run, J2 n >o ^ s = N wnere s — 2™, because, in 
the end, all sites are wetted. For RG, one expects 
P s ~ s 2 ~ T ~ s -5 ' 91 , so P s is a slowly decreasing func- 
tion of s, until s approaches the size of the system, at 
which point the "infinite" clusters contribute. 

In Fig. [IT] we plot the average of the normalized dis- 
tribution, (P s )/N, as a function of s, for different L. For 
the RG case, the data show expected decrease with s, 
except for a large accumulation in the last two bins be- 
cause of the finite size effects. On the other hand, for the 
PR case, the P s /N seem to be increasing, except pos- 
sibly for a small region in the largest system, and the 
accumulation in the large bins occurs over a much wider 
range. 

In the lower plots of Fig. [TT] we show the slopes be- 
tween pairs of points (taking the logarithm of (P s )/N 
first). The data of the slopes for the RG model for large 
L is seen to be consistent with the theoretical prediction 
2 — t = — 5/91. For the PR model, for smaller s and L, 
the slopes are positive, consistent with the observation 
of [2] who found r = 1.9(1). Of course, for a normaliz- 
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FIG. 9. (Color online) Horizontal wrapping probability. Expanded plots are given in the lower panels, which show the precise 
crossing of the curves for different L; note that the horizontal scale for the PR model is more expanded than that of the RG 
model. 
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FIG. 10. (Color online) Estimates of p c vs. L -0 ' 75 (RG) and L _1 (PR). In each case, we have estimates determined from the 
average of p at which one-way wrapping occurs (top), the average value of p at which horizontal wrapping first occurs (middle), 
and the average value of p at which either horizontal or vertical crossing first occurs (bottom). Linear fits to the data are 
shown in the figure, where y represents the p c estimates and x represents 1/L. 



able size distribution, it is necessary that r > 2, at least 
asymptotically. We indeed find that the slope (barely) 
goes below zero (corresponding to r > 2) for a range 
of s for the largest system; however, it is unclear from 
these data whether the slope truly approaches a consis- 
tent value or whether it contains for example logarithmic 
terms. Simulations on larger systems should help to an- 
swer this question. 

For the corrections to scaling for the critical size dis- 



tribution, one expects 

P s ~ s 2 - T {A + Bs- n ...) . (3) 

The data for RG are consistent with Vl w 0.75 as found 
previously [551 • If we fit the data of the PR model to 
we find B is negative, Q « 0.3, and r w 2.025. The 
latter value is consistent with /3/u = 0.05 through the 
scaling relation t = 2 + j3/(i/D), assuming D = 2. The 
hyper-scaling relation j3 /v = d— D would imply that D sa 
1.95, and the scaling is also consistent with this value 
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FIG. 11. (Color online) {P s )/N vs. InL (upper plots), and logarithmic slopes between pairs of points (lower plots), for L - 
2048, and 8192 (peaks going from left to right). 



512, 



of D. Thus, there is evidence that the size distribution 
becomes power-law and that scaling is satisfied for the 
situations in which v ^ 1 and 7/1/ 7^ 2. 



IV. CONCLUSIONS 

We have found the critical bond fraction p c for the PR 
model on the square lattice to high accuracy by a number 
of methods. The two best criteria to determine p c (in 
terms of convergence with L) are the average value of p 
at which horizontal wrapping first occurs (Fig. 10 1, and 
the crossing point of the horizontal wrapping probability 
(Fig. [9]). (One could just as well use vertical wrapping, 
or the sum of the two |30j . as a criterion.) Combining 
our measurements, we conclude 



Pc = 0.526565 ± 0.000005 



(4) 



where the error bars represent a combination of statis- 
tical error and also the variation among results based 
upon different criteria. This is consistent with the value 
0.5266(2) given by Radicchi and Fortunato [j"4"] . 

The striking qualitative different between the explosive 
and regular percolation is highlighted by the non-zero 
limiting behaviors of M'^pj/N (Fig. [I]) and x/iV (Fig. 
[5]). These results suggest a discontinuity at the transi- 
tion point, in contrast to RG (regular percolation), where 



the corresponding quantities are continuous. (Note that 
for regular percolation on a hierarchical small-world net- 
work, however, the transition can also be discontinuous 

EH). 

The cluster size distribution of the PR model shows 
quite different behavior than the ER model, with possible 
power-law behavior for very large systems with strong 
finite-size effects. 

The wrapping probabilities proved useful for locat- 
ing the transition point and, perhaps surprisingly, be- 
have qualitatively quite similar to the RG model. The 
horizontal wrapping probability shows a very well- 
defined crossing point, just as found for the RG case. 
Its value, n^- 1 = 0.5160, is quite close to (but not iden- 
tical with) the value for standard percolation, IL^ — 
0.52105829 . . .[53]. This result recalls the recent findings 
of various kinetic systems that evolve to mimic random 
percolation [3"21 13"3"] . 

On the other hand, the value of the one-way wrapping 
probability, E[W for the PR model (see Fig. is quite a 
bit below the RG percolation value, 0.351642855 . . ., and 
it is hard to find its asymptotic value precisely Evidently, 
because of the more compact geometry of the PR giant 
cluster, wrapping one way is more difficult than in the 
RG case. 

Finally, for the scaling, we have found some contradic- 
tory results: M2, (s max ), IlW, and the size distribution 
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give P/v = 0.06(1), 7/1/ = 1.90(1), and v = 1.04(1), and 
t = 2.025(10), implying D = 1.94(1), where number in 
parentheses represents our estimated errors in the last 
digit (s), while some of the other results (such as the be- 
havior of M' 2 and £) are more consistent with v = 1 and 
j/p = 2. Perhaps this is indicative that the normal two- 
parameter scaling does not hold for this model because of 
the first-order transition, or that logarithmic corrections 
come into play. 

Note added: While this paper was in revision, a 
preprint appeared which argues that the explosive perco- 
lation transition in the case of the PR rule on the random 



graph is continuous [34] . Those arguments however do 
not apply to the regular square lattice studied here. 
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